O objetivo desse notebook é aplicar a Linguagem R em Ciência de Dados na resolução de problemas de classificação.
A base de dados iris fornece medidas em centímetros das variáveis largura e comprimento da sépala e largura e comprimento da pétala, respectivamente, para 50 flores, de cada uma das 3 espécies de íris.
data('iris')
head(iris)
Podemos usar o scatterplot para entender o comportamento das classes
plot(iris[,1:4], col=iris$Species)
É muito importante tratar valores ausêntes. A função is.na e is.nan nos diz se um elemento é NA ou NaN. Combinada com a função any, podemos descobrir se algum elemento da base é NA.
print(any(is.na(iris)))
## [1] FALSE
print(any(apply(iris, 2, is.nan)))
## [1] FALSE
Devemos normalizar os dados, uma vez que existem classificadores que são suscetíveis a escala dos atributos.
iris = data.frame(scale(iris[,1:4]), Species=iris$Species)
summary(iris)
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Min. :-1.86378 Min. :-2.4258 Min. :-1.5623 Min. :-1.4422
## 1st Qu.:-0.89767 1st Qu.:-0.5904 1st Qu.:-1.2225 1st Qu.:-1.1799
## Median :-0.05233 Median :-0.1315 Median : 0.3354 Median : 0.1321
## Mean : 0.00000 Mean : 0.0000 Mean : 0.0000 Mean : 0.0000
## 3rd Qu.: 0.67225 3rd Qu.: 0.5567 3rd Qu.: 0.7602 3rd Qu.: 0.7880
## Max. : 2.48370 Max. : 3.0805 Max. : 1.7799 Max. : 1.7064
## Species
## setosa :50
## versicolor:50
## virginica :50
##
##
##
Depois de normalizar a base, precisamos particionar o conjunto de treinamento e teste. Para essa e outras bases vamos usar o holdout com 70% para treinamento e 30% para teste.
set.seed(123)
idx = sample(nrow(iris), 0.7*nrow(iris), replace=FALSE)
tran = iris[idx,]
test = iris[-idx,]
Podemos conferir se ocorreu contaminação.
intersect(rownames(tran), rownames(test))
## character(0)
Agora podemos induzir alguns modelos de Aprendizado de Máquina. Vamos começar pela Árvore de Decisão.
require('rpart')
## Loading required package: rpart
model = rpart(Species~., iris)
No caso da Árvore de Decisão, é possível plotar o modelo. Para isso vamos usar a biblioteca rattle.
library('rattle')
## Loading required package: tibble
## Loading required package: bitops
## Rattle: A free graphical interface for data science with R.
## Version 5.4.0 Copyright (c) 2006-2020 Togaware Pty Ltd.
## Type 'rattle()' to shake, rattle, and roll your data.
fancyRpartPlot(model, caption=NULL)
Também podemos predizer a classe das amostras de teste.
pred = predict(model, test, type='class')
E plotar o scatterplot para os dados de treinamento e teste.
par(mfrow=c(1,2))
plot(tran[,3:4], col=tran$Species, pch=as.numeric(tran$Species))
plot(test[,3:4], col=pred, pch=as.numeric(pred))
Podemos avaliar de forma mais objetiva o desempenho do modelo.
conf = table(test$Species, pred)
print(conf)
## pred
## setosa versicolor virginica
## setosa 14 0 0
## versicolor 0 18 0
## virginica 0 1 12
Inclusive, calcular a taxa de acerto.
acc = sum(diag(conf))/sum(conf)
print(acc)
## [1] 0.9777778
Outras meedidas de desempenho poderiam ser utilizadas. Recomendo olhar os pacotes ‘caret’ e ‘mlr’ para mais informaões..
#F1-sconre
f1 <- function(conf) {
prc = diag(conf)/colSums(conf)
rec = diag(conf)/rowSums(conf)
aux = (2 * prc * rec) / (prc + rec)
return(mean(aux))
}
#AUC
require('ROCR')
## Loading required package: ROCR
auc <- function(class, pred) {
aux = prediction(as.vector(pred), as.vector(class))
aux = unlist(performance(aux, "auc")@y.values)
return(aux)
}
Agora podemos induzir outros modelos de Aprendizado de Máquina. Para isso precisamos carregar algumas bibliotecas.
require('kknn')
## Loading required package: kknn
require('e1071')
## Loading required package: e1071
require('randomForest')
## Loading required package: randomForest
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:rattle':
##
## importance
Gerando os modelos de Aprendizado de Máquina.
model1 = kknn(Species~., tran, test, k=3)
model2 = svm(Species~., tran, kernel="radial")
model3 = randomForest(Species~., tran)
Calculando a predição.
pred1 = model1$fitted.values
pred2 = predict(model2, test, type='class')
pred3 = predict(model3, test)
Calculando o desempenho baseado na acurácia.
acuracia <- function(orig, pred) {
acc = table(orig, pred)
sum(diag(acc))/sum(acc)
}
cat('A acurácia do 3NN é', acuracia(test$Species, pred1))
## A acurácia do 3NN é 0.9333333
cat('A acurácia do SVM é', acuracia(test$Species, pred2))
## A acurácia do SVM é 0.9777778
cat('A acurácia do RF é', acuracia(test$Species, pred3))
## A acurácia do RF é 0.9777778
No entanto, o holdout não é a melhor alternativa. A validação cruzada estratificada é mais interessante. Podemos ver na imagem um exemplo de seu funcionamento.
k-fold cross validation (Adaptado da Wikepedia)
Podemos implementar :)
kfolds <- function(y, folds) {
names(y) = 1:length(y)
index = lapply(1:nlevels(y), function(i) {
rep(1:folds, length.out=length(y[y == levels(y)[i]]))
})
index = unlist(index)
folds = lapply(1:folds, function(i) {
as.numeric(names(y[index == i]))
})
return(folds)
}
E agora podemos usar.
idx = kfolds(iris$Species, 10)
Antes vamos conferir se ocorreu contaminação.
intersect(idx[[1]], idx[[2]])
## numeric(0)
E agora, de fato, podemos usar.
acc = lapply(idx, function(i){
model = randomForest(Species~., iris[-i,])
pred = predict(model, iris[i,])
acuracia(iris[i,]$Species, pred)
})
cat(mean(unlist(acc)), '+-', sd(unlist(acc)))
## 0.9533333 +- 0.03220306
Um experimento interessante é a tarefa de classificar dígitos utilizando o algoritmo de Redes Neurais Artificiais. Para isso precisamos relizar os experimentos em duas fases: (1) treinar a rede para reconhecer os dígitos e (2) utilizá-la para classificar novos dígitos. Para isso precisamos carregar as blbiotecas ‘tensorflow’ e ‘keras’.
require('tensorflow')
## Loading required package: tensorflow
require('keras')
## Loading required package: keras
A base de dados de dígitos se chama MNIST e pode ser carregada com a função dataset_mnist. Ela contém 60000 amostras dos dígitos 0 até 9. Cada dígito é uma imagem de 28 x 28 pixels.
mnist = dataset_mnist()
A função show_digit recebe um dígito na forma de uma matriz 28 x 28 e imprime o dígito na forma de imagem.
show_digit <- function(numero) {
image(numero[1:28,28:1])
}
Podemos utilizar a função show_digit para plotar uma amostra de todos os dígitos presentes na base de dados.
par(mfrow=c(2,5))
plot = lapply(0:9, function(i) {
aux = mnist$train$x[mnist$train$y == i,,][1,,]
show_digit(t(aux))
})
Normalizando os dados de treinamento e teste.
mnist$train$x = mnist$train$x/255
mnist$test$x = mnist$test$x/255
Construindo uma Rede Neural Artificial usando Keras. Essa rede é composta de uma camanda de entrada flatten (achatada de 28 x 28 neurônios), uma camada oculta dense (128 neurônios com função de ativação linear) e uma camada de saída dense (10 neurônios com função de ativação exponencial normalizada).
modelo = keras_model_sequential() %>%
layer_flatten(input_shape = c(28, 28)) %>%
layer_dense(units = 128, activation = "relu") %>%
layer_dense(10, activation = "softmax")
Informações sobre a Rede Neural Artificial.
summary(modelo)
## Model: "sequential"
## ________________________________________________________________________________
## Layer (type) Output Shape Param #
## ================================================================================
## flatten (Flatten) (None, 784) 0
## ________________________________________________________________________________
## dense_1 (Dense) (None, 128) 100480
## ________________________________________________________________________________
## dense (Dense) (None, 10) 1290
## ================================================================================
## Total params: 101,770
## Trainable params: 101,770
## Non-trainable params: 0
## ________________________________________________________________________________
Podemos representar graficamente:
library(deepviz)
library(magrittr)
modelo %>% plot_model()
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
Ainda precisamos definir o algoritmo de treinamento, a função de avaliação e a medida de desempenho que vamos utilizar. Podemos fazer isso com a função compile.
modelo %>%
compile(
loss = "sparse_categorical_crossentropy",
optimizer = "adam",
metrics = "accuracy"
)
Agora podemos treinar a rede.
historico = modelo %>%
fit(
x = mnist$train$x, y = mnist$train$y,
epochs = 10,
validation_split = 0.3,
verbose = 2
)
É interessante medir o erro e o desempenho da rede ao longo das épocas no conjunto de treinamento e validação.
plot(historico)
## `geom_smooth()` using formula 'y ~ x'
Agora que sabemos que a rede aprendeu a reconhecer dígitos, podemos verificar seu desempenho no conjunto de teste. Para isso precisamos aplicar essas amostras nunca antes vistas na entrada da rede e obter os valores de saída.
predicao = predict(modelo, mnist$test$x)
head(predicao)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 9.558057e-09 3.779262e-11 1.319685e-08 9.883004e-05 3.488893e-14
## [2,] 4.910057e-12 9.836884e-07 9.999982e-01 6.844302e-07 8.030510e-19
## [3,] 2.539843e-07 9.997995e-01 4.072102e-06 4.887717e-07 7.108847e-07
## [4,] 9.999865e-01 1.239132e-11 3.391276e-06 5.745591e-10 2.942153e-11
## [5,] 2.212367e-09 4.571540e-14 8.091424e-08 9.510420e-11 9.996395e-01
## [6,] 6.279258e-09 9.999075e-01 1.821460e-09 4.694487e-09 3.349559e-08
## [,6] [,7] [,8] [,9] [,10]
## [1,] 1.054294e-08 2.430504e-17 9.999007e-01 3.572829e-07 7.341557e-08
## [2,] 7.384209e-08 1.229754e-10 3.823800e-17 4.582389e-08 8.587397e-18
## [3,] 1.153944e-05 4.967835e-06 1.623248e-04 1.592033e-05 2.178347e-07
## [4,] 2.806128e-06 3.090005e-06 4.101109e-06 1.618421e-11 9.907203e-08
## [5,] 1.237713e-09 1.869342e-09 9.177735e-06 1.342951e-07 3.510211e-04
## [6,] 2.322132e-09 3.439437e-09 9.235813e-05 7.486649e-08 1.311218e-08
Cada neurônio de saída da rede irá retornar um valor. Aquele neurônio com maior valor deve indicar a classe com maior proababilidade.
predicao = apply(predicao, 1, which.max) -1
predicao[1:6]
## [1] 7 2 1 0 4 1
Calculando o desempenho da rede no conjunto de teste.
acc = acuracia(mnist$test$y, predicao)
cat('A acurácia do modelo no conjunto de teste eh:', acc)
## A acurácia do modelo no conjunto de teste eh: 0.976
Agora podemos usar os modelos de Árvore de Decisão, kNN, SVM e Random Forest para comparar com o desempenho com a Rede Neural Artificial. Primeiro vamos construir uma base tabular com a mnist.
tran = data.frame(mnist$train$x, y=as.factor(mnist$train$y))
test = data.frame(mnist$test$x, y=as.factor(mnist$test$y))
Agora podemos induzir outros modelos.
model = rpart(y~., tran)
pred = predict(model, test, type='class')
cat('A acurácia da ANN é:', acc)
## A acurácia da ANN é: 0.976
cat('A acurácia da DT é', acuracia(test$y, pred))
## A acurácia da DT é 0.6196